zsim = 189*ones(Nsim,Tsim);
ssim = 189*ones(Nsim,Tsim);
gsim = 189*ones(Nsim,Tsim);
nsim = 189*ones(Nsim,Tsim);
bsim = 189*ones(Nsim,Tsim+1); ibsim = 189*ones(Nsim,Tsim+1);
rsim = 189*ones(Nsim,Tsim);
fsim = 189*ones(Nsim,Tsim+1);
dsim = 189*ones(Nsim,Tsim);


%-initialize
itshow = 0;
for inn = 1:Nsim
    idef = 0;  % start with market access

    if itshow <= 0 
        disp(['sim ', num2str(inn),'/',num2str(Nsim)])
        itshow = 500;
    end
    itshow = itshow - 1;

    for tt = 1:Tsim
        if tt == 1
            ibtp = ib0; igtp = ig0; istp = is0; iztp = iz0; 
            bsim(inn,tt) = bvec(ibtp);
            fsim(inn,tt) = fsim0;
        else
             run draw_random_state.m
        end  
    
        zsim(inn,tt) = zvec(iztp);
        ssim(inn,tt) = istp;
        gsim(inn,tt) = gvec(igtp);
        
        if idef == 0 % market access this period
            
            dsim(inn,tt) = dpol(ibtp,igtp,istp,iztp);
    
            if dsim(inn,tt) == 0 % no default
                ibnp = bopt_loc(ibtp,igtp,istp,iztp);
                bsim(inn,tt+1)  = bpol(ibtp,igtp,istp,iztp);
                ibsim(inn,tt+1) = bopt_loc(ibtp,igtp,istp,iztp);
                nsim(inn,tt)    = nsched(ibnp,ibtp,igtp,istp);
                rsim(inn,tt)    = rhosched(ibnp,igtp,istp);
                
                fsim(inn,tt+1)  = (1.0D0-dlta)*(fsim(inn,tt)/gsim(inn,tt)) + (nsim(inn,tt)/gsim(inn,tt));
                idefnp          = 0;
            
            elseif dsim(inn,tt) == 1  % default this period            
                bsim(inn,tt+1)  = bsim(inn,tt)/gsim(inn,tt);
                ibsim(inn,tt+1) = -1;
                nsim(inn,tt)    = 0;
                rsim(inn,tt)    = -8;
                idefnp          = 1;
            else
                disp(['something wrong: dpol = ', num2str(dsim(inn,tt))])
                messi
            end
        
        elseif idef == 1 % no market access this period
            % should add re-entry here     
            dsim(inn,tt) = 1;
            idefnp       = 1;
        end
        
    
        %--udpate state
        iglp = igtp;
        islp = istp;
        idef = idefnp;
        ibtp = ibsim(inn,tt+1);
    end
end

%% statistics

IRF_rr = 189*ones(Tsim,1); IRF_rr_sd = 189*ones(Tsim,1); 
IRF_nn = 189*ones(Tsim,1); IRF_nn_sd = 189*ones(Tsim,1); 
IRF_bb = 189*ones(Tsim,1); IRF_bb_sd = 189*ones(Tsim,1); 
IRF_ff = 189*ones(Tsim,1); IRF_ff_sd = 189*ones(Tsim,1); 
for tt = 1:Tsim
    dd = dsim(:,tt); xx = 1-dd;
    rr = rsim(:,tt); nn = nsim(:,tt); bb = bsim(:,tt); ff = fsim(:,tt);

    rrx = rr(xx==1); nnx = nn(xx==1); bbx = bb(xx==1); ffx = ff(xx==1);

    IRF_rr(tt,1) = mean(rrx); IRF_rr_sd(tt,1) = std(rrx);
    IRF_nn(tt,1) = mean(nnx); IRF_nn_sd(tt,1) = std(nnx);
    IRF_bb(tt,1) = mean(bbx); IRF_bb_sd(tt,1) = std(bbx);
    IRF_ff(tt,1) = mean(ffx); IRF_ff_sd(tt,1) = std(ffx);

end